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Abstract 



j^ I We present the results of large scale computer simulations in which we in- 

c/3 ' vestigate the static and dynamic properties of sodium disilicate and sodium 

trisilicate melts. We study in detail the static properties of these systems, 
namely the coordination numbers, the temperature dependence of the Q'"-* 
species and the static structure factor, and compare them with experiments. 
We show that the structure is described by a partially destroyed tetrahe- 



O ' dral Si04 network and the homogeneously distributed sodium atoms which 

are surrounded on average by 16 silicon and other sodium atoms as nearest 
neighbors. We compare the diffusion of the ions in the sodium silicate systems 
^ i with that in pure silica and show that it is much slower in the latter. The 

C^^ \ sodium diffusion is characterized by an activated hopping through the Si-0 

matrix which is frozen with respect to the movement of the sodium atoms. 

^D ■ We identify the elementary diffusion steps for the sodium and the oxygen 

diffusion and find that in the case of sodium they are related to the breaking 

Q> , of a Na-Na bond and in the case of oxygen to that of a Si-0 bond. From 

the self part of the van Hove correlation function we recognize that at least 
two successive diffusion steps of a sodium atom are spatially highly correlated 
with each other. With the same quantity we show that at low temperatures 
also the oxygen diffusion is characterized by activated hopping events. 



I. INTRODUCTION 



;_i ' Silicate melts and glasses are an important class of materials in very different fields, 

e.g. in geosciences (since silicates are geologically the most relevant materials) and in tech- 
nology (windows, containers, and optical fibers). From a physical point of view it is a very 
challenging task to understand the properties of these materials on a microscopic level, and 
in the last twenty years many studies on different systems have shown that molecular dy- 
namics computer simulations are a very appropriate tool for this purpose (Angell et ai, 1981, 
Balucani and Zoppi, 1994, Kob, 1999). The main advantage of such simulations is that they 
give access to the whole microscopic information in form of the particle trajectories which 
of course cannot be determined in real experiments. 



In pure silica (Si02) the structure is that of a disordered tetrahedral network in which 
Si04 tetrahedra are connected via the oxygens at the edges. In a recent simulation (Horbach 
and Kob, 1999a) we have studied in detail the statics and dynamics of pure silica. In the 
present paper we investigate the statics and dynamics of melts in which the network modifier 
sodium is added to a silica melt. In particular two sodium silicate systems are discussed, 
namely sodium disilicate (Na2Si205) and sodium trisilicate (Na2Si307) to which will be 
referred to NS2 and NS3, respectively. 

In recent years several authors found, by using the potential proposed by Vessal et al. 
(1989) that, e.g., NS2 is characterized by a microsegregation in which the sodium atoms 
form clusters of a few atoms between bridged Si04 units (Vessal et a/., 1992, Smith et a/., 
1995, Cormack and Cao, 1997). This result is somewhat surprising because in experiments 
one observes a phase separation between Na20 and Si02 at lower temperatures and more 
likely in sodium silicates with a higher Si02 concentration (Mazurin et al., 1970, Haller et 
al., 1974). In order to see whether microsegregation is also reproduced with a different model 
from the one of Vessal et al. we have performed our simulations with a different potential 
(discussed below). 

Up to now the dynamics of systems like alkali silicates was investigated by studying 
only the dynamics of the alkali atoms which, at lower temperatures, is much faster than 
that of the silicon and oxygen atoms. In this paper we therefore present the results for the 
diffusion dynamics of all components in NS2 and NS3 as well as the microscopic mechanism 
of diffusion in these systems. 

II. MODEL AND DETAILS OF THE SIMULATIONS 

In a classical molecular dynamics computer simulation one solves numerically Newton's 
equations of motion for a many particle system. If quantum mechanical effects can be ne- 
glected such simulations are able to give in principle a realistic description of any molecular 
system. The determining factor of how well the properties of a real material are reproduced 
by a MD simulation is the potential with which the interaction between the atoms is de- 
scribed. The model potential we use to compute the interaction between the ions in sodium 
silicates is the one proposed by Kramer et al. (1991) which is a generalization of the so called 
BKS potential (van Beest et al., 1990) for pure silica. It has the following functional form: 

0(r) = ^^ + A„^exp(-i?«^r)-% a,/5 G [Si, Na, O]. (1) 

Here r is the distance between an ion of type a and an ion of type (3. The values of the 
parameters A^p, B^p and C^/? can be found in the original publication. The potential (p has 
been optimized by Kramer et al. for zeolites, i.e. for systems that have Al ions in addition 
to Si, Na and O. In that paper the authors used for silicon and oxygen the partial charges 
gsi = 2.4 and go = —1-2, respectively, whereas sodium was assigned its real ion charge 
gNa = 1-0. With this choice charge neutrality is not fulfilled in sodium silicate systems. To 
overcome this problem we introduced for the sodium ions a distance dependent charge g(r) 
instead of gNa, 

,. r0.6(l + ln[C(re-r)' + ll) r<r, ,^. 

O.D r > Tr 



which means that for r > r^ charge neutrahty is vahd {q{r) = 0.6 for r > Tc). Note that 
q{r) is continuous at Tc. We have fixed the parameters r^ and C such that the experimental 
mass density of NS2 and the static structure factor from a neutron scattering experiment 
are reproduced well. From this fitting we have obtained the values Tc = 4.9 A and C = 
0.0926 A~^. With this choice the charge g(r) crosses smoothly over from g(r) = 1.0 at 1.7 
A to q{r) = 0.6 for r > Tc. 

The simulations have been done at constant volume with the density of the system fixed 
to 2.37 g/cm^. The systems consist of 8064 and 8016 ions for NS2 and NS3, respectively. 
The reason for using such a relatively large system size is that, like in Si02 (Horbach et al, 
1996, Horbach et al, 1999b), strong finite size effects are present in the dynamics of smaller 
systems which have to be avoided (Horbach and Kob, 1999c). The equations of motion were 
integrated with the velocity form of the Verlet algorithm and the Coulombic contributions 
to the potential and the forces were calculated via Ewald summation. The time step of the 
integration was 1.6 fs. The temperatures investigated are 4000 K, 3400 K, 3000 K, 2750 K, 
2500 K, 2300 K, 2100 K, and in addition for NS2 also 1900 K. The temperature of the system 
was controlled by coupling it to a stochastic heat bath, i.e. by substituting periodically 
the velocities of the particles with the ones from a Maxwell-Boltzmann distribution with 
the correct temperature. After the system was equilibrated at the target temperature, we 
continued the run in the microcanonical ensemble, i.e. the heat bath was switched off. In 
order to improve the statistics we have done two independent runs at each temperature. 
The production runs have been up to 7.5 ns real time which corresponds to 4.5 million time 
steps. In the temperature range under investigation the pressure decreases monotonically 
from 3.8 GPa to 1.8 GPa in the case of NS3 and from 7.3 GPa to 4.1 GPa in the case of NS2. 
In order to make a comparison of the static structure factor for NS2 from our simulation 
with one from a neutron scattering experiment (see below) we have also determined the 
structures of the glass at T = 300 K. The glass state was produced by cooling the system 
from equilibrated configurations at T = 1900 K with a cooling rate of 1.16 ■ 10^^ K/s. The 
pressure of the system at T = 300 K is 0.96 GPa. 

In the following we compare the properties of NS2 and NS3 with those of pure silica 
which we have investigated in recent simulation. The details of the latter can be found in 
Horbach and Kob (1999a). We only mention here that these simulations were done also at 
the density 2.37g/cm in the temperature range 6100 K > T > 2750 K for a system of 8016 
particles. The two production runs at the lowest temperature were over about 20 ns real 
time, and the pressure at this temperature is 0.9 GPa. 

III. RESULTS 

A. Structural properties 

In order to investigate the local environment of an atom in a disordered structure, espe- 
cially at high temperatures, it is useful to calculate the coordination number Zapir) which 
gives the number of atoms of type j3 e [Si, Na, O] surrounding an atom of type a G [Si, Na, O] 
within a distance r' < r. Note that Za/^ir) is essentially the integral from r' = to r' = r 
over the function 47rr'^5f„^ (r') , where Qapif') denotes the pair correlation function for a/3 
correlations (Balucani and Zoppi, 1994). In Fig. ^a we show -2si-o('") and -Zo-Si('") for Si02 



at T = 2750 K and for NS3 and NS2 at T = 2100 K. In both quantities we observe a strong 
increase for distances from 1.5 A to 1.7 A at which the functions reach a plateau which per- 
sists up to about 3.0 A. For distances r > 3.0 A there is no such step hke behavior because of 
the disorder. The vertical lines in Fig. [^a at fsi_o = 1-61 A and r^^^ = 2.35 A correspond 
to the first maximum and the first minimum in gsio{f), respectively. In zsi-ol^) the plateau 
is essentially at z = 4 in Si02 as well as in the sodium silicate systems which means that 
most of the silicon atoms, i.e. more than 99 %, are surrounded by four oxygen atoms thus 
forming a tetrahedron. In the case of Si02 there is also a plateau essentially at ^o-Si = 2 
in the same r interval. Therefore, in this system most of the the oxygen atoms are bridging 
oxygens between two tetrahedra, and thus the system seems to form a perfect disordered 
tetrahedral network even at the relatively high temperature T = 2750 K. Indeed we found 
that at this temperature less than one percent of the silicon and oxygen atoms are defects 
(Horbach and Kob, 1999a). For NS3 and NS2 a plateau is formed at -^o-Si = 1-71 and 
2;o-si = 1-6, respectively. This is due to the fact that only a part of the oxygen atoms are 
bridging oxygens, namely 68.5 % in the case of NS2 and 71.3 % in the case of NS3. Most of 
the other oxygen atoms form dangling bonds with one silicon neighbor (28.4 % in NS2 and 
22.7 % in NS3 at r^\n°) or they are not nearest neighbors of silicon atoms (3.1 % in NS2 
and 6.0 % in NS3 at rj[^ ). The local environment of the sodium atoms is characterized 
by 2Na-o(^)) ^Na-si('")5 and 2;Na-Na('^) in Fig. |l|b. In these quantities there is no step like 
behavior also at small distances. At r^fn"°, which is approximately at r = 3.0 A for NS3 
and NS2, the apparent coordination number z = 4.2 is too small in comparison with X-ray 
experiments from which one would expect a value between 5 and 6 (see Brown et al, 1995, 
and references therein). Also fNa-o = 2.2 A is too small in comparison with the values found 
in X-ray experiments which are between 2.3 A and 2.6 A. We recognize in Fig. |l|b that in 
the case of NS2 the functions ZNa-Na('") and 2Na-si('") are relatively close to each other for 
r < r^tn^^'^r^tn^' ~ 5.1 A. So at r = r^f^"^^ and r = r^t~^' the apparent coordination 
numbers are 2;Na-Na = 7.8 and ^Na-si = 8.6, respectively, i.e., are quite similar. In contrast 
to that, the coordination numbers in NS3 at r = r^fn"^*^ and r = r^fn"^' are ;ZNa-Na = 5.8 
and 2;Na-Si('^) = 9.8, respectively, which means that there is a substitution effect in NS3 
such that more silicon atoms are on average in the neighborhood of an sodium atom than 
in NS2 because there are less sodium atoms in NS3. Therefore, in NS2 as well as in NS3 
every sodium atom has on average about 16 silicon and sodium atoms in its neighborhood. 
The structure of NS2 and NS3 can be studied in more detail by looking at the so called 
Q*-"^ species which can be determined experimentally by NMR (Stebbins, 1995) and by 
Raman spectroscopy (McMillan and Wolf, 1995). Q^"-* is defined as the fraction of Si04 
tetrahedra with n bridging oxygens in the system. In Fig. |T]a we have seen that these 
structural elements are very well defined also at temperatures as high as T = 2100 K. From 
Fig. 1^ we recognize that at T = 2100 K there is essentially a ternary distribution of Q^'^\ 
Q^'^\ and Q^^^ both in NS2 and in NS3. At higher temperatures there is also a significant 
contribution of silicon defects, i.e. three- and five-fold coordinated silicon atoms (denoted 
as rest in the figure), and of Q^^^ species. For T < 3200 K the curves for Q^'^^ are well 
described by Arrhenius laws /(T) oc exp{Eji^/T) (bold solid lines in Fig. §) with activation 
energies E^ = 5441 K and Ea = 6431 K for NS2 and NS3, respectively. If we extrapolate 
these Arrhenius laws to low temperatures we recognize from Fig. ^ that the Q*-^-* species 
essentially disappear slightly above the experimental glass transition temperatures Tg exp for 



NS2 and NS3, which are at 740 K and 760 K (Knoche et ai, 1994), respectively. Therefore, 
if we would be able to equilibrate our systems down to Tg exp we would expect a binary 
distribution of Q^^^ and Q*-^-* species in the glass state of NS2 and NS3. Also included in 
Fig. ^ is the Q^"^ species distribution for NS2 at T = 300 K which we have calculated from 
the configurations which we have produced by cooling down the system from T = 1900 K 
to T = 300 K with a cooling rate of 10^^ K/s. We recognize from the data that the Q^"^ 
species distribution essentially coincides at T = 300 K with the one at T = 2100 K, which is 
due to the fact that the liquid structure at the latter temperature has just been frozen in. If 
we extrapolate the data from 2100 K to lower temperatures allowing further relaxation we 
would expect about 45 % Q^^^ and 55 % Q^'') in the case of NS3 and about 40 % Q^^^ and 
60 % Q^'^^ in the case of NS2 at Tg^exp- In contrast to these values one finds in NMR and 
Raman experiments 60 % Q^^^ and 40 % Q^^^ in the case of NS3 and 92 % Q^^) ^nd 8 % Q^^^ 
in the case of NS2 (Stebbins, 1988, Mysen and Frantz, 1992, Sprenger et al, 1992, Knoche, 
1993). Thus our simple model is not able to describe the Q'-"-' species reliably. The reason 
for this is probably that our model gives a slightly wrong coordination function -^Na-ol'") as 
we have mentioned before. 

So far we have looked at the local environment of the atoms. In order to investigate the 
structure on a larger length scale useful quantities are the partial static structure factors, 

Sapiq) = ^ im (exp(iq ■ (r^ - r„))) , (3) 

^^ 1=1 m=l 

depending on the magnitude of the wave-vector q. The factor /^/j is equal to 0.5 for a 7^ /3 
and equal to 1.0 for a = p. As examples Figs, ^ja and ^d show Sooil) and 5'NaNa(?) 
at the temperatures T = 4000 K and T = 2100 K for NS2 and NS3, respectively. The 
peak at 2.8 A~^ in Soo{q) for NS2 and NS3 corresponds to the length scale 27r/2.8 A^^ = 
2.24 A which is approximately the period of the oscillations in gooii^)- In the same way 
the peak at 2.1 A~^ in S^aNa{c[) is due to the period of oscillations in 5'NaNa('")- Moreover, 
the temperature dependence we observe for the peaks at 2.8 A~^ and 2.1 A~^ is relatively 
weak. At T = 2100 K a peak is visible at 1.7 A^^ in Soo{q) for NS2 and NS3, which is 
the so called first sharp diffraction peak (FSDP). This feature arises from the tetrahedral 
network structure since the length scale which corresponds to it, i.e. 27r/1.7 A = 3.7 A, 
is approximately the spatial extent of two connected Si04 tetrahedra. Only a shoulder can 
be identified in Soo{q) for T = 4000 K at the position of the FSDP because the network 
structure is less pronounced than at T = 2100 K. In agreement with the aforementioned 
interpretation of the FSDP this feature is absent in the Na-Na correlations. We note that 
this holds also for the Na-0 and Si-Na correlations. We see therefore that the structure 
in NS2 is characterized by a partially destroyed tetrahedral network, and that on the other 
hand there are the sodium atoms which are homogeneously distributed in the system having 
on average 16 sodium and silicon neighbors. We would therefore expect that there is a 
characteristic length scale of regions where the network is destroyed and where the sodium 
atoms are located. This explains the peak at g = 0.95 A"^ in the Safsiq) corresponding to 
a length scale 27r/0.95 A = 6.6 A which is approximately two times the mean distance of 
nearest Na-Na or Na-Si neighbors. We mention that this peak is also present in the Si-Si, 
Si-0, Si-Na, and Na-0 correlations. We emphasize that no evidence for a microsegregation 



is found in the partial structure factors both for NS2 and for NS3, at least at the density of 
this study 2.37 g/cm . 

In order to make a comparison of the static structure factor for NS2 measured in a 
neutron scattering experiment by Misawa et al. (1980) with that from our model we have 
to determine S^^^{q) which is calculated by weighting the partial structure factors from the 
simulation with the experimental coherent neutron scattering lengths ha {a G [Si,Na, O]): 

S''''''{q) = ^^^Y.bahpSaM- (4) 

Z^o ^^a'^a af3 

The values for 6„ are 0.4149 ■ lO'^^ ^^^ o.363 ■ lO'^^ ^^ g^^^ 0.5803 ■ 10"^^ cm for silicon, 
sodium and oxygen, respectively. They are taken from Susman et al. (1991) for silicon and 
oxygen and from Bacon (1972) for sodium. Fig. ^ shows S^'^^{q) from the simulation and 
the experiment at T = 300 K. We see that the overall agreement between simulation and 
experiment is good. For q > 2.3 A~^ the largest discrepancy is at the peak located at 
q = 2.8 A~^ where the simulation underestimates the experiment by approximately 15% in 
amplitude. Very well reproduced is the peak at g = 1.7 A^^ The peak at g = 0.95 A^^ is 
not present in the experimental data which might be due to the fact that this feature is less 
pronounced in real systems. 



B. Dynamical properties 

One of the simplest quantities to study the dynamics of liquids is the self diffusion 
constant Da for a tagged particle of type a, which can be calculated from the mean squared 
displacement (r^(t)) via the Einstein relation. 

The different Da for NS2 and NS3 are shown in Fig. ^ as a function of the inverse temper- 
ature. Also included are the diffusion constants for pure silica from our recent simulation 
(Horbach and Kob, 1999a). For the latter system Dsi and Dq show at low temperatures 
the expected Arrhenius dependence, i.e. Da ex: exp [EA/{kBT)]. The corresponding acti- 
vation energies Ea are in very good agreement with the experimental values by Brebec et 
al. (1980) and Mikkelsen (1984) (given in the figure). We recognize from Fig. ^ that the 
dynamics in the sodium silicate melts is much faster than in Si02 even at the relatively high 
temperature T = 2750 K (10^/T = 3.64 K~^) for which D^i and Dq are about two orders 
of magnitude larger in NS2 and NS3 than in Si02. Furthermore, we see that the sodium 
diffusion constants can be fitted in the whole temperature range very well with Arrhenius 
laws. The corresponding activation energies are Ea = 0.93 eV for NS2 and Ea = 1.26 eV 
for NS3. These activation energies are about 20 % to 30 % higher than those obtained from 
electrical conductivity measurements (Greaves and Ngai, 1995, and references therein). But 
this discrepancy might at least be partly due to the fact that our simulations are done at 
relatively high pressures. With decreasing temperature DNa decouples more and more from 
-Do and Dsi in NS2 and NS3. Therefore, at least at low temperatures, the motion of the 
oxygen and silicon atoms is frozen with respect to the timescale of motion of the sodium 
atoms. 



Experimentally the diffusion constants for oxygen and silicon have been measured by 
Poe et al. (1997) for the system Na2Si409 at high temperatures and high pressures. Besides 
a few other combinations of temperature and pressure they determined Z)si and Dq at 
T = 2800 K and p = 10 GPa where they found the values Dgi = 1-22 • 10^ cmVs and Dq = 
1.53 ■ 10^ cm^/s. In order to make a comparison with these values we have done a simulation 
of Na2Si409 for a system of 7680 particles at the density 2.9 g/cm^. At this density we have 
found at T = 2800 K and p = 10.13 GPa the diffusion constants Dsi = 0.78 ■ 10^ cmVs, 
Dq = 1.17 ■ 10^ cm^/s, and DNa = 2.35 ■ 10^ cm^/s. Thus D^i and Dq from our simulation 
underestimate the experimental data by less than 40% which is within the error bars of 
the experiment. Therefore it seems that our model gives a quite realistic description of the 
diffusion in sodium silicate melts. 

In order to give insight into the microscopic mechanism which is responsible for the 
diffusion of the different species in NS2 and NS3 we discuss now the time dependence of 
the probability Pa/sit) that a bond between an atom of type a and an atom of type (3 
is present at time t when it was present at time zero. For this we define two atoms as 
bonded if their distance from each other is less than the location of the first minimum rmin 
in the corresponding partial pair correlation function gapif)- In the following we restrict our 
discussion to NS2 because the conclusions which are drawn below also hold for NS3. From 
the functions Qapir) for NS2 we find for rmin the values 3.6 A, 5.0 A, 2.35 A, 5.0 A, 3.1 A, and 
3.15 A for the Si-Si, Si-Na, Si-0, Na-Na, Na-0, and 0-0 correlations. As an example Fig. |^ 
shows PNa-Na for the different temperatures. First of all we recognize from this figure that a 
plateau is formed on an intermediate time scale which becomes more and more pronounced 
the lower the temperature is. This plateau is due to the fact that at r = rmin the amplitude 
of fifNa-Nal^) IS equal 0.56 and not zero. Thus there are some sodium atoms which vibrate 
between the first and the second neighbor shell leading to the first fast decay of PNa-Na to the 
plateau value. In the long time regime of PNa-Na, where it decays from the plateau to zero, 
its shape seems to be independent of temperature. This is also true for the functions Papif) 
for the other correlations. For this reason it makes sense to define the lifetime Tap of a bond 
between two atoms of type a and (3 as the time at which Papit) has decayed to 1/e. Indeed, if 
the function PNa-Na for the different temperatures is plotted versus the scaled time t/rNa-Na 
(inset of Fig. ^) one obtains one master curve in the long time regime t/rNa-Na > 0.1. 
This master curve cannot be described by an exponential function but is well described by a 
Kohlrausch- Williams-Watts (KWW) function, P(t) ex exp(— (t/rNa-Na)^) with an exponent 
f3 = 0.54, which is shown in Fig. ^ in which this function is fitted to the curve for T = 1900 K. 

The lifetimes r^/j can now be correlated with the diffusion constants by plotting different 
products Tap ■ D^ versus temperature, which is done in Fig. [^. The product TNa-Na ■ -^Na is 
essentially constant over the whole temperature range whereas TNa-o ■ -^Na increases with 
decreasing temperature. This means that the elementary diffusion step for the sodium 
diffusion is related to the breaking of an Na-Na bond and not to that of an Na-0 bond, 
although the nearest neighbor distance is smaller for Na-0 (rNa-o = 2.2 A) than for Na-Na 
(^Na-Na = 3.3 A). We meutiou that an Arrhenius law holds also for TNa-o for which we have 
found the activation energy E^ = 1.14 eV. The latter can be interpreted as an effective Na- 
O binding energy in NS2. In NS3 TNa-o can be fitted with an Arrhenius law for T < 3000 K, 
the binding energy in this case is 1.64 eV. Also constant is the product Tsi_o ■ -Do which 
shows that the oxygen diffusion is related to the breaking of Si-0 bonds. In contrast to that 



rsi_o ■ -Dsi is only constant at high temperatures. For temperatures T < 3000 K this product 
decreases with decreasing temperature. Concerning the oxygen and sihcon diffusion we have 
found recently that the same conclusions hold also for pure silica (Horbach and Kob, 1999a). 

We have seen that the temperature dependence of Dq and Dsi for Si02 changes strongly 
with the addition of sodium atoms and, moreover, diffusion becomes much faster. In Fig. |^ 
we compare the behavior of the quantity Psi-o(^) for NS2, NS3 and Si02 at T = 2750 K. 
We recognize that the shape of the curves seems to be the same for the three systems. The 
only difference lies in the time scale which is, as expected from the behavior of the diffusion 
constants, about two orders of magnitudes larger for silica than the sodium silicate systems. 
That the shape of the curves is indeed the same for the three systems is demonstrated in 
the inset of Fig. ^ in which we have plotted the same data as before versus the scaled time 
t/rsi-o- We see that the three curves fall nicely onto one master curve which can be fitted 
by a KWW law with exponent j3 = 0.9. This is an astonishing result for Psi-o(^) since the 
local environment of the oxygen atoms is very different in the sodium silicate systems from 
that in silica. 

Further insight on how diffusion takes place in sodium silicates can be gained from the 
self part of the van Hove correlation function which is defined by (Balucani and Zoppi, 1994) 

1 A^a 

G:{r,t) = —J2{5ir-\r,{t)-rM\))) «G{Si,Na,0} . (6) 

Thus 47rr^G"(r, t) is the probability to find a particle a distance r away from the place it 
was at t = 0. In Figs. ^ and |pD we show this probability for different times at T = 2100 
K for sodium and oxygen, respectively. Note that we have chosen in both cases a linear- 
logarithmic plot. At t = 0.6 ps the sodium function exhibits a single peak with a shoulder 
around r > 2.8 A. This shoulder becomes more pronounced as time goes on and at t = 6.7 
ps there is a second peak located at a distance r which is equal to the average distance 
^Na-Na = 3.3 A bctwecu two nearest sodium neighbors (marked with a vertical line in 
Fig. I^a). The first peak is still located at the same position as it was at t = 0.6 ps while 
its amplitude has decreased. Thus we can conclude from this that the sodium atoms do 
not diffuse in a continuous way but discontinuously in time by hopping on average over 
the distance fNa-Na- This interpretation was first given for similar features in a soft-sphere 
system by Roux et al. (1989). At t = 45.7 ps even a third peak has developed at 2fNa-Na 
while the first two peaks remain at the same position. This means that many sodium 
atoms have performed now a second diffusion step. We see from this that two successive 
elementary diffusion steps, each corresponding to a breaking of a Na-Na bond, are spatially 
highly correlated with each other. At t = 164.5 ps the function has lost its three peak 
structure but the first peak is still observable with a significant amplitude. 

Whereas many of the sodium atoms have performed two elementary diffusion steps at 
t = 45.7 ps the amplitude of the first peak for oxygen (Fig. |b) decreases only from about 
1.5 to 0.8 in the time interval 0.6 ps < t < 45.7 ps. In this time interval most of the oxygen 
atoms sit in the cage which is formed by the neighboring atoms and only rattle around in 
this cage. Nevertheless, in this time window a shoulder becomes more and more pronounced 
around the mean distance between two nearest oxygen neighbors fo-o = 2.61 A. This means 
that also the oxygen diffusion takes place by activated hopping events. We have found the 
same behavior for oxygen at low temperatures in pure silica (Horbach and Kob, 1999a). 

8 



IV. SUMMARY 

By using a simple pair potential we have performed large scale molecular dynamics sim- 
ulations in order to investigate the dynamic properties of the two sodium silicate melts 
Na2Si205 (NS2) and Na2Si307 (NS3). The structure of these two systems can be character- 
ized by a partially destroyed tetrahedral Si04 network, and a homogeneous distribution of 
sodium atoms which have on average about 16 sodium and silicon atoms as nearest neigh- 
bors. The regions in between the network structure introduce a new length scale which is 
about two times the distance of nearest Na-Na or Na-Si neighbors. This leads to a peak in 
the static structure factor at g = 0.95 A^^. Furthermore, we have found no evidence in the 
static structure factor that a microsegregation of Na20 complexes takes place. A comparison 
of experimental data to the one of our computer simulation shows that the latter gives a fair 
description of the static properties of NS2 and NS3. We have explicitly demonstrated this 
by showing that the static structure factor S'^'^"^{q) from a neutron scattering experiment for 
NS2 (Misawa et ai, 1980) is reproduced quite well by our simulation. In contrast to our 
simulation this experiment exhibits no peak at g = 0.95 A~^ which might be due to the fact 
that this feature is less pronounced in real systems. Our simulations give not a very good 
description for the Q*^"-* species distribution which is perhaps due to the fact that our model 
underestimates the coordination number of nearest Na-0 neighbors. 

Nevertheless the oxygen and silicon diffusion constants for Na2Si409 are in very good 
agreement with experimental data by Poe et al. (1997). In comparing the diffusion constants 
in NS2 and NS3 with those in silica we have recognized that the dynamics becomes much 
faster with the addition of the sodium atoms. Moreover, in the sodium silicates the diffusion 
constant for sodium decouples for decreasing temperature more and more from those of 
oxygen and silicon, such that at low temperatures the dynamics of the silicon and oxygen 
atoms is frozen in with respect to the movement of the sodium atoms. The sodium diffusion 
constants for NS2 and NS3 exhibit over the whole temperature range we investigate an 
Arrhenius behavior with activation energies around 1 eV. 

We have shown that our simulation is able to give insight into the microscopic mecha- 
nism of diffusion in NS2 and NS3. From the time dependent bond probability Papit) we 
determined an average lifetime r^^ of bonds between an atom of type a and an atom of type 
p. By correlating the lifetimes Tajs with the diffusion constants, we show that the elemen- 
tary diffusion step in the case of sodium is the breaking of a Na-Na bond and in the case of 
oxygen that of a Si-0 bond. By studying the self part of the van Hove correlation function 
we have demonstrated that at low temperatures the diffusion for oxygen and sodium takes 
place by activated hopping events over the mean distance fo-o = 2.6 A and fNa-Na = 3.3 A, 
respectively. Moreover, we have found for the case of sodium that at least two successive 
diffusion steps are spatially highly correlated with each other. 
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FIG. 1. Coordination numbers 2^/3 as a function of distance r at the temperature T = 2750 K 
for Si02 and at T = 2100 K for NS2 and NS3. a) zsi-o{r) and zo-si{r), b) ZNa-o(^)) •2Na-Si('")) 
and ZNa-Na(^)- For the explanation of the vertical lines see text. 
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FIG. 2. Temperature dependence of the Q^"'' species for NS2 (closed symbols) and NS3 (open 
symbols). The bold lines are fits with Arrhenius laws with the activation energies Ea = 5441 K 
for NS2 and Ea = 6431 K for NS3. 
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FIG. 3. Partial static structure factors Soo{q) (thin lines) and <S'NaNa('?) (bold lines) at the 
temperatures T = 4000 K (dashed lines) and T = 2100 K (sohd lines), a) NS2, b) NS3. 



13 




0.0 2.0 4.0 6.0 8.0 10.0 

° —1 



FIG. 4. Comparison of the static structure factor S^^^{q) from our simulation (solid line) with 
the experimental data of Misawa et al. (1980) (dashed line). 
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FIG. 5. The diffusion constants for Si02, NS2, and NS3. The bold straight lines are fits with 
Arrhenius laws. The experimental values of the activation energies for silica are taken from Brebec 
et al. (1980) for silicon and from Mikkelsen (1984) for oxygen. 
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FIG. 6. Time dependence of PNa-Naj the probability that a bond between two sodium atoms 
which exists at time zero is also present at time t, for all temperatures investigated. Inset: Plot of 
the same data versus the scaled time t/iTMa-Na as a function of temperature. 
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FIG. 7. The products r^/? • Da show whether or not the diffusion constant Da is correlated with 
the lifetime of a bond r^^. 
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FIG. 8. Psi-o(i) for NS2, NS3, and Si02 at T = 2750 K. Inset: Plot of the same data versus 
the scaled time t/rsi_o- 
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FIG. 9. Space- and time-dependence of the self part of the van Hove correlation function for 
NS2 at T = 2100 K in a linear-logarithmic plot, a) for sodium atoms, b) for oxygen atoms. The 
vertical lines correspond to fNa-Na = 3.3 A and 2fNa-Na in a) and fo-o = 2.61 A in b). 



16 



